#####################################################
#####################################################
###Sensitivity Analyses (GLM, LM and ZINNB Models)###
#####################################################
#####################################################
library(doBy)
library(ggplot2)
library(mgcv)
library(dplyr)
library(interplot)
library(lfe)
library(survival)
library(MASS)
library(foreign)
library(stargazer)
library(pscl)
library(sandwich)
library(lmtest)
###Set working directory
setwd("~/OneDrive - Indiana University/FromGoogle/DeforestationDisease/CodeFinal_12_5_24/")
full.2003.2019<- read.csv("afrogriddiseasejaxa20032019.csv")
summary(full.2003.2019$year)
###Create deforestation indicators
#Contemp
full.2003.2019$deforestation_shock2 <- ifelse(full.2003.2019$jaxa_forest_y_n==1, full.2003.2019$NDVI_mean_anom, 0)
summary(full.2003.2019$deforestation_shock2)

#####################
#####################
###Linear Models. ###
#####################
#####################

##########
###GLMs###
##########

###Socioeconomic model
pois_glm_sep1 <- glm(confirm ~ deforestation_shock2+
                       log.nl + logged.pop, data=full.2003.2019, family="poisson")
summary(pois_glm_sep1)
#Clustered SEs
pois_glm_sep1.cl <- coeftest(pois_glm_sep1, vcov. = cluster.vcov(pois_glm_sep1, full.2003.2019$gid))

#Full model
pois_glm_sep2 <- glm(confirm ~ deforestation_shock2+
                       log.nl + logged.pop +
                       logged_state + logged_nonstate +
                       p_anom+spei+t_anom+f_month2+f_month3+
                       f_month4+f_month5+f_month6+f_month7+f_month8+
                       f_month9 +f_month10+f_month11+
                       f_month12, data=full.2003.2019, family="poisson")
summary(pois_glm_sep2)
#Clustered SEs
pois_glm_sep2.cl <- coeftest(pois_glm_sep2, vcov. = cluster.vcov(pois_glm_sep2, full.2003.2019$gid))

##Optimized model
pois_glm_sep3 <- glm(confirm ~ deforestation_shock2+
                       log.nl + logged.pop + logged_nonstate + 
                       t_anom + f_month5 + f_month6 + f_month7 + f_month9, 
                     data=full.2003.2019, family="poisson")
summary(pois_glm_sep3)
#Clustered SEs
pois_glm_sep3.cl <- coeftest(pois_glm_sep3, vcov. = cluster.vcov(pois_glm_sep3, full.2003.2019$gid))
pois_glm_sep3.cl

##Optimized RE model
pois_glm_sep4 <- glmer(confirm ~ deforestation_shock2+
                         log.nl + logged.pop + logged_nonstate +
                         t_anom + f_month5 + f_month6 + f_month7 + f_month9+
                         (1|gid)+(1|month),
                       data=s.dat2, family="poisson")
summary(pois_glm_sep4)

##Adjust AIC for missingness
#Subset NA data
s.dat <- na.omit(full.2003.2019[,c("confirm", "deforestation_shock2","log.nl","logged.pop","logged_state","logged_nonstate","p_anom","spei","t_anom","f_month2","f_month3","f_month4","f_month5","f_month6","f_month7","f_month8","f_month9","f_month10","f_month11","f_month12","c_gid","gid","month")])
#Socioeconomic model
pois_glm_sep1s <- glm(confirm ~ deforestation_shock2+
                        log.nl + logged.pop, data=s.dat, family="poisson")
summary(pois_glm_sep1s)
#Clustered SEs
pois_glm_sep1s.cl <- coeftest(pois_glm_sep1s, vcov. = cluster.vcov(pois_glm_sep1s, s.dat$gid))

#Optimized model
pois_glm_sep3s <- glm(confirm~deforestation_shock2+
                        log.nl + logged.pop + logged_nonstate + 
                        t_anom + f_month5 + f_month6 + f_month7 + f_month9, 
                      data=s.dat, family="poisson")
summary(pois_glm_sep3s)
#Clustered SEs
pois_glm_sep3s.cl <- coeftest(pois_glm_sep3s, vcov. = cluster.vcov(pois_glm_sep3s, s.dat$gid))

#Optimized RE model
pois_glm_sep4s <- glmer(confirm ~ deforestation_shock2+
                          log.nl + logged.pop + logged_nonstate +
                          t_anom + f_month5 + f_month6 + f_month7 + f_month9+
                          (1|gid)+(1|month),
                        data=s.dat, family="poisson")
summary(pois_glm_sep4s)


#Export to LaTex
stargazer(pois_glm_sep1.cl,pois_glm_sep2.cl,pois_glm_sep3.cl)
stargazer(pois_glm_sep1,pois_glm_sep2,pois_glm_sep3)

#Compare AICs
pois_glm_sep1$aic
pois_glm_sep2$aic
pois_glm_sep3$aic
AIC(pois_glm_sep4)
pois_glm_sep1s$aic
pois_glm_sep3s$aic
AIC(pois_glm_sep4s)



##########################
###Overdispersion tests###
##########################
###Socioeconomic model
#Subset NA rm data
m.dat <- na.omit(full.2003.2019[,c("confirm", "deforestation_shock2","log.nl","logged.pop","c_gid")])
#1. Create fitted values
Phats1<-fitted.values(pois_mod_sep1)
#2. Caluclate u_i
Uhats1<-((m.dat$confirm-Phats1)^2 - m.dat$confirm) / (Phats1 * sqrt(2))
#3. Estimate a linear regression to see if we can reject the null
summary(lm(Uhats1~Phats1))
#So, we do not reject null, suggesting overdispersion does not exists

###Full model
#1. Create fitted values
Phats2<-fitted.values(pois_mod_sep2)
#2. Caluclate u_i
Uhats2<-((s.dat$confirm-Phats2)^2 - s.dat$confirm) / (Phats2 * sqrt(2))
#3. Estimate a linear regression to see if we can reject the null
summary(lm(Uhats2~Phats2))
#So, we do not reject null, suggesting overdispersion does not exists

###Optimized model
#Subset NA rm data
opt.dat<- na.omit(full.2003.2019[,c("confirm", "deforestation_shock2","log.nl","logged.pop","logged_nonstate","t_anom","f_month5","f_month6","f_month7","f_month9")])

#1. Create fitted values
Phats3<-fitted.values(pois_mod_sep3)
#2. Caluclate u_i
Uhats3<-((opt.dat$confirm-Phats3)^2 - opt.dat$confirm) / (Phats3 * sqrt(2))
#3. Estimate a linear regression to see if we can reject the null
summary(lm(Uhats3~Phats3))
#So, we do not reject null, suggesting overdispersion does not exists`



#########
###AMs###
#########
####Socioeconomic model
am1 <- gam(confirm ~ s(deforestation_shock2, bs = "cr")+
                  log.nl + logged.pop+c_gid, data=full.2003.2019)
summary(am1)

####Full
am2 <- gam(confirm ~ s(deforestation_shock2, bs = "cr")+
                  log.nl + logged.pop +
                  logged_state + logged_nonstate +
                  p_anom+spei+t_anom+f_month2+f_month3+
                  f_month4+f_month5+f_month6+f_month7+f_month8+
                  f_month9 +f_month10+f_month11+
                  f_month12+c_gid, data=full.2003.2019)
summary(am2)

####Optimized model
am3 <- gam(confirm ~ s(deforestation_shock2, bs = "cr")+
             log.nl + logged.pop +
             logged_nonstate +
             t_anom+f_month5+f_month6+f_month7+
             f_month9, data=full.2003.2019)
summary(am3)

####Optimized RE model
am4 <- gam(confirm ~ s(deforestation_shock2, bs = "cr")+
             log.nl + logged.pop +
             logged_nonstate +
             t_anom+f_month5+f_month6+f_month7+
             f_month9+
             s(gid, bs = "re")+s(month, bs = "re"), data=full.2003.2019)
summary(am4)


###Adjust AIC for full model N
###Socioeconomic model
am1s <- gam(confirm ~ s(deforestation_shock2, bs = "cr")+
                   log.nl + logged.pop+c_gid, data=s.dat)
summary(am1s)
##Optimized model
am3s <- gam(confirm ~ s(deforestation_shock2, bs = "cr")+
             log.nl + logged.pop +
             logged_nonstate +
             t_anom+f_month5+f_month6+f_month7+
             f_month9, data=s.dat)
summary(am3s)
####Optimized RE model
am4s <- gam(confirm ~ s(deforestation_shock2, bs = "cr")+
             log.nl + logged.pop +
             logged_nonstate +
             t_anom+f_month5+f_month6+f_month7+
             f_month9+
             s(gid, bs = "re")+s(month, bs = "re"), data=s.dat)
summary(am4s)

##Compare AICs
am1$aic
am2$aic
am3$aic
AIC(am4)
#Adjusted for N
am1s$aic
am3s$aic
AIC(am4s)

#########
###LMs###
#########
###Socioeconomic model
lm1 <- lm(confirm ~ deforestation_shock2+
                  log.nl + logged.pop, data=full.2003.2019)
summary(lm1)
#Clustered GID estimates
lm1.cl <- coeftest(lm1, vcov. = vcovCL(lm1, full.2003.2019$gid))

#Full model
lm2 <- lm(confirm ~ deforestation_shock2+
                  log.nl + logged.pop +
                  logged_state + logged_nonstate +
                  p_anom+spei+t_anom+f_month2+f_month3+
                  f_month4+f_month5+f_month6+f_month7+f_month8+
                  f_month9 +f_month10+f_month11+
                  f_month12, data=full.2003.2019)
summary(lm2)
#with(summary(lm2), 1 - deviance/null.deviance)
#Clustered GID estimates
lm2.cl <- coeftest(lm2, vcov. = vcovCL(lm2, full.2003.2019$gid))

####Optimized model
lm3 <- lm(confirm ~ deforestation_shock2+
            log.nl + logged.pop +
            logged_nonstate +
            t_anom+f_month5+f_month6+f_month7+
            f_month9, data=full.2003.2019)
summary(lm3)
#with(summary(lm2), 1 - deviance/null.deviance)
#Clustered GID estimates
lm3.cl <- coeftest(lm3, vcov. = vcovCL(lm3, full.2003.2019$gid))

####Optimized RE model
lm4 <- lmer(confirm ~ deforestation_shock2+
              log.nl + logged.pop + logged_nonstate +
              t_anom + f_month5 + f_month6 + f_month7 + f_month9+
              (1|gid)+(1|month),
            data=full.2003.2019)
summary(lm4)

##Adjust AIC for missingness
##Socioeconomic model
lm1s <- lm(confirm ~ deforestation_shock2+
                   log.nl + logged.pop, data=s.dat)
summary(lm1s)
#Clustered GID estimates
lm1s.cl <- coeftest(lm1s, vcov. = vcovCL(lm1s, s.dat$gid))

##Optimizied model
lm3s <- lm(confirm ~ deforestation_shock2+
             log.nl + logged.pop +
             logged_nonstate +
             t_anom+f_month5+f_month6+f_month7+
             f_month9, data=s.dat)
summary(lm3s)
#Clustered GID estimates
lm3s.cl <- coeftest(lm3s, vcov. = vcovCL(lm3s, s.dat$gid))

####Optimized RE model
lm4s <- lmer(confirm ~ deforestation_shock2+
              log.nl + logged.pop + logged_nonstate +
              t_anom + f_month5 + f_month6 + f_month7 + f_month9+
              (1|gid)+(1|month),
            data=s.dat)
summary(lm4s)


#Export to LaTex
stargazer(lm1.cl,lm2.cl, lm3.cl)
stargazer(lm1,lm2,lm3)

#AIC
AIC(lm1)
AIC(lm2)
AIC(lm3)
AIC(lm4)
#Adjust AIC
AIC(lm1s)
AIC(lm3s)
AIC(lm4s)

##################################
##################################
###Additional Robustness Models###
##################################
##################################

#################
###Country FEs###
#################
###Socioeconomic model
cfe_sep1 <- felm(confirm ~ deforestation_shock2+
                   log.nl + logged.pop|
                   COWCODE|0|gid, data=full.2003.2019)
summary(cfe_sep1)
#Full model
cfe_sep2 <- felm(confirm ~ deforestation_shock2+
                   log.nl + logged.pop +
                   logged_state + logged_nonstate +
                   p_anom+spei+t_anom+f_month2+f_month3+
                   f_month4+f_month5+f_month6+f_month7+f_month8+
                   f_month9 +f_month10+f_month11+
                   f_month12|
                   COWCODE|0|gid, data=full.2003.2019)
summary(cfe_sep2)
#Optimized model
cfe_sep3 <- felm(confirm ~ deforestation_shock2+
                   log.nl + logged.pop +logged_nonstate +
                   t_anom+f_month5+f_month6+f_month7+
                   f_month9|
                   COWCODE|0|gid, data=full.2003.2019)
summary(cfe_sep3)

###Run models adjusted for N
#NA omit
s.dat2 <- na.omit(full.2003.2019[,c("confirm", "deforestation_shock2","log.nl","logged.pop","logged_state","logged_nonstate","p_anom","spei","t_anom","f_month2","f_month3","f_month4","f_month5","f_month6","f_month7","f_month8","f_month9","f_month10","f_month11","f_month12","COWCODE","gid")])
##Soc.
cfe_sep1s <- felm(confirm ~ deforestation_shock2+
                    log.nl + logged.pop|
                    COWCODE|0|gid, data=s.dat2)
##Opt.
cfe_sep3s <- felm(confirm ~ deforestation_shock2+
                    log.nl + logged.pop +logged_nonstate +
                    t_anom+f_month5+f_month6+f_month7+
                    f_month9|
                    COWCODE|0|gid, data=s.dat2)


##Export to Latex
stargazer(cfe_sep1,cfe_sep2,cfe_sep3)

###AIC
##CFE models
AIC(cfe_sep1)
AIC(cfe_sep2)
AIC(cfe_sep3)
#adjusted
AIC(cfe_sep1s)
AIC(cfe_sep3s)

##########################
##########################
###Zero Inflated Models###
##########################
##########################

###Socioeconomic model
zip_glm_sep1 <- zeroinfl(confirm ~ deforestation_shock2+
                           log.nl + logged.pop|
                           log.nl + logged.pop+logged_state + logged_nonstate +
                           p_anom+spei+t_anom+factor(month), data=full.2003.2019, family="poisson")
summary(zip_glm_sep1)

#Full model
zip_glm_sep2 <- zeroinfl(confirm ~ deforestation_shock2+
                           log.nl + logged.pop +
                           logged_state + logged_nonstate +
                           p_anom+spei+t_anom|
                           log.nl + logged.pop+ logged_state + logged_nonstate +
                           p_anom+spei+t_anom+factor(month), data=full.2003.2019, family="poisson")
summary(zip_glm_sep2)

#Optimized model
zip_glm_sep3 <- zeroinfl(confirm ~ deforestation_shock2+
                           log.nl + logged.pop +logged_nonstate +
                           t_anom+f_month5+f_month6+f_month7+
                           f_month9|
                           log.nl + logged.pop+logged_state + logged_nonstate +
                           p_anom+spei+t_anom+factor(month), data=full.2003.2019, family="poisson")
summary(zip_glm_sep3)


##Export to LaTex
stargazer(zip_glm_sep1,zip_glm_sep2,zip_glm_sep3)
stargazer(zip_glm_sep1,zip_glm_sep2,zip_glm_sep3, zero.component = T)
AIC(zip_glm_sep1,zip_glm_sep2,zip_glm_sep3)





